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As the universe evolves, it becomes more inhomogeneous due to gravitational clumping. We 
attempt to find a function that characterizes this behavior and increases monotonically as inhomo- 
geneity increases. We choose S — In f2 as the candidate "gravitational entropy" function, where f2 is 
the phase-space volume below the Hamiltonian H of the system under consideration. We perform 
a direct calculation of Q for transverse electromagnetic waves and gravitational waves, radiation 
and density perturbations in an expanding FLRW universe. These calculations are carried out in 
the linear regime under the assumption that the phases of the oscillators comprising the system are 
random. Entropy is thus attributed to the lack of knowledge of the exact field configuration. The 
time dependence of H leads to a time-dependent f2. We find that Q, and hence InQ behaves as 
required. We also carry out calculations for Bianchi IX cosmological models and find that, even 
in this homogeneous case, the function can be interpreted sensibly. We compare our results with 
Penrose's C 2 hypothesis. Because S is defined to resemble the fundamental statistical mechanics def- 
inition of entropy, we are able to recover the entropy in a variety of familiar circumstances including, 
evidently, black-hole entropy. The results point to the utility of the relativistic ADM Hamiltonian 
formalism in establishing a connection between general relativity and statistical mechanics, although 
fully nonlinear calculations will need to be performed to remove any doubt. 



I. INTRODUCTION 



It has been recognized for some time that gravity behaves in an"antithermodynamic" fashion. Whereas ordinary 
thermodynamic systems, a gas for example, tend to become more homogenous with time, gravitating systems tend 
to become more inhomogeneous with time. The anomalous behavior can be viewed as a manifestation of the long- 
range nature of the gravitational force, which tends to cause the components of a gravitating system to clump. If 
we associate an increase in homogeneity with an increase in entropy for thermodynamic systems, then for gravitating 
systems an increase in entropy will imply an increase in inhomogeneity. The "gravitational arrow of time" points in 
the direction of increasing inhomogeneity. 

There have, apparently, been only a few attempts in the literature to characterize the gravitational arrow of time. 
The most well-known is the suggestion of Penrose [Q that "gravitational entropy" should be measured by C 2 , the 
square of the Weyl tensor. Penrose hoped that the Weyl tensor would provide a measure of inhomogeneity and increase 
monotonically in time. This proposal met with some degree of success with a slight modification Q-H]. However, this 
"entropy" function is not well defined for all spacctimcs; for example conformally flat or vacuum models. Furthermore 
Bonnor p] has found an example in which the gravitational arrow points in the opposite sense when compared to 
the flow of radiation from a collapsing fluid, throwing doubt on the entire proposal. There have been several other 
efforts to define the entropy of the gravitational field from various standpoints (Smolin Q , Hu and Kandrup JtJ and 
Brandenberger et al. Q), but none appear to have established an explicit connection to the Hamiltonian formulation 
of gravity, and none has addressed the arrow-of-time question. 

In this paper we attack the problem of gravitational entropy by a direct approach. The goal is to find a function 
that behaves like entropy, i.e. that increases monotonically as a gravitating system becomes more inhomogeneous. 
We therefore choose a function that resembles entropy as much as possible: 

5 = In ft (1) 

Here, S is gravitational entropy and £1 is the volume of phase space for the system. (Unless stated otherwise, 
throughout the paper we use units in which h — c — k ~ G = 1). 

For this choice we have reverted to the fundamental statistical mechanics definition of entropy. However, although 
we will refer to particle models, it is absolutely crucial to realize our goal is to characterize the phase space and 
entropy of the field itself, not of systems of particles. 

There are several advantages and disadvantages to the above definition for gravitational entropy. For thermodynamic 
systems, a direct evaluation of the phase space is extremely difficult, if not impossible. Instead, one chooses the simpler 
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path of evaluating the partition function, Z = J^i e , from which the entropy is readily derived as S = k(\n Z+(3E) , 
where (3 = l/(kT) and E is the mean energy. 

Here, however, we encounter the first conceptual difficulty in carrying over the procedure to relativity: To evaluate 
the partition function requires knowing the temperature of the system. In general relativity we usually deal with 
dynamical, not thermodynamic, systems, and a temperature is not well-defined. A macroscopic pendulum executing 
simple harmonic motion, for example, constitutes a dynamical, not a thermodynamic system. Of course, one could 
assign an effective temperature kT ~ mv 2 where v is the pendulum's velocity, but the system is nevertheless not in 
thermal equilibrium and so the concept of a partition function is not obviously useful. 

However, a pendulum's motion does define a trajectory in phase space and it can be calculated without recourse 
to temperature. This is one of the two main reasons for reverting to the statistical mechanics expression for entropy 
(|l|). The other, anticipating the application to cosmology, is that the phase space approach is intimately connected 
with Hamiltonian dynamics, and a Hamiltonian formalism of relativity (the ADM formalism ||) already exists. 

Now, normally, one is not interested in the absolute value of the entropy, but in the change of entropy with time. 
If the above pendulum were given a higher energy, it would include a larger phase space and the logarithm of the 
area between the two paths would formally resemble an entropy change. However, neglecting dissipative forces, the 
pendulum does not change its trajectory. Furthermore, the concept of entropy implies a loss of information, i.e. that 
we do not know the pendulum's location within this band. This requires that the system be ergodic, which is not 
true in the case of the pendulum. 

On the other hand, if one imagines a system with N independent oscillators and assumes that their trajectories are 
uncorrelated, in other words, that their phases are random, then one should be able to compute an entropy via ([!]). 

Thus our approach is basically simple: we calculate the phase space of dynamical systems in relativity, assuming 
that the trajectories of the components are independent and that consequently each region of phase space is occupied 
with equal probability. We then derive an entropy via (|]) . (The number of degrees of freedom need not be large as 
long as the system is chaotic, as in the case of Bianchi IX cosmologies.) 

Usually, one computes the entropy as the logarithm of a volume of phase space constrained by an energy E. 
Our approach, when applied to cosmological models, forces us to substitute the Hamltonian H for E. This is the 
most natural and conservative extention of the usual definition but it should be emphasized that H does not always 
correspond to the energy. In most of the systems we consider, H will be time-dependent, resulting in an entropy 
change. 

The main advantage of the entire method is its conceptual clarity. The main disadvantage of the procedure is that it 
is technically cumbersome. However, we have found a number of systems for which the computation is tractable in the 
classical, perturbative limit. In these limits the entropy function does appear to increase or decrease monotonically 
when appropriate, and by suitable identification of parameters we recover the entropy familiar from a variety of 
circumstances including, evidently, black holes. In this way the universality of the phase-space concept is established. 
Extensions to the nonpertubative and quantum limits need to be carried out. 

In §§[n] — IV of the paper we summarize some preliminary calculations necessary for what follows. In §^ we apply 
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IX. In §|X| we compare our method with Penrose's C 2 suggestion. In §pOTwe argue that our entropy is indeed the 
Bekenstein-Hawking entropy under appropriate circumstances, and in § KII| we discuss further applications of our 
procedure. 



our meth od to the electromagnetic field. In §VI and § VII we treat gravitationa l wa ves and density perturbations. 



In § VIII| w e give a more formal mathematical basis for the results of §|VI| and §|VII. Section |X| concerns Bianchi 



II. TWO PARTICLE MODEL 



We now present a simple model to illustrate the basic approach. This model describes the Newtonian gravitational 
interaction between two particles, and although it is an extremely idealized particle model, it does highlight several 
important aspects of our treatment that will remain unchanged in the more complex field models. 

Consider two particles, each of mass m free to slide on a 1-dimensional frictionless track of length L with hard 
"bumpers" set at the two ends. The Hamiltonian for this system is 

H-lL + lL + V^), (2) 



with 



Gm 2 

V = - . „ , ■ (3) 
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The factor ro softens the potential and is introduced to avoid singularities. The system thus represents two objects 
that can pass through each other, such as colliding galaxies. In this case the Hamiltonian H — Eq, is the total energy, 
which remains constant. 

The "air-track" model is a closed, isolated system and the available phase space can be computed. For two particles, 
however, the system cannot be regarded as ergodic and hence an entropy is not really defined. When generalized to 
N 1 particles and three dimensions the system may be regarded as a microcanonical ensemble and the particles 
can be assumed to be in random motion and equally likely to be found in any region of phase space. In that case 
an entropy would be well-defined. Unfortunately, the complexity of a 2N dimensional phase space prohibits analytic 
solutions; hence we demonstrate the basic results with just two particles then argue a connection to more general 
iV-body systems. 

For the two-particle case, we can write the phase space below energy Eq as 



ft(E < E ) 



L/2 <.x 2 „ 
-L/2 •' x 2 m , 



-^2m(E-V) f\f2m(E-V)-p\ 

dx2 I dpi / dp2 



-y/2m(E-V) J - y f2m(E-V)-p 2 1 



(4) 



Note this procedure is similar to evaluating the volume of a 4-sphere, although the exact topology and hence volume 
will depend on the form of V. 

Generally one defines the acessible region of phase space as a shell between Eq and Eq + AE; thus £1 = FI(Eq < E < 
(Eq + AE)). It is, however, easier to evaluate the full volume fl(E < Eq), a procedure that we will follow throughout 
the paper. In the limit of a large number of degrees of freedom, the two results are identical, since most of the volume 
of an TV-object resides infinitesimally near the surface. 

The first important step in the phase-space procedure is to find the limits of integration, which are not always 
obvious. Due to the quadratic form of the momenta in the Hamiltonian (^|), the momentum integrals give the volume 
of a 2-sphere and the limits are set simply by requiring the pf to be positive definite but constrained by the total 
energy of the system, as in (0). After evaluating the momentum integrals we have 



r L/2 rx 2maI 

Q(E < E ) = 27rm / dx x / (E — V)dx 2 - 

-L/2 J%2 min 



(5) 



The lower and upper limits on x 2 are set by restricting our attention to bound systems, such that E < and by 
requiring E — V > 0, which leads to 



X2 



Xt, ± 



G 2 m 4 



E 2 'o- 

The entire volume of phase space can then be evaluated analytically and is found to be 



n(E < E Q ) = 4irm 3 GL 



sinh 





(6) 



(7) 



where Vq = —Gm 2 /ro is the minimum potential. 

Note several aspects of this result. For a fixed form of the potential V, there are only two ways to change the phase 
space: one must change either Eq or Vq. As expected, for a larger Eq, the particles are free to roam around in a larger 
region of phase space and thus fl increases. Were E to decrease due to dissipation, the particles would be confined 
to a smaller volume. This is one example of gravitational clumping. Now, to change Vq, one must change ro, the 
softening parameter. Decreasing ro makes the potential deeper and vice-versa. Thus, imposing a finite ro has the 
effect of excluding a certain region of phase space compared to the usual gravitational potential, in which V — > — oo 
as r — > 0. The dependence of Q on Eq and ro is shown more explicitly in Fig. ^ where we plot phase space trajectories 
for the two particle system, assuming one particle to be stationary with zero momentum. 

The parameter Eq merely determines whether the system is overall bounded. For the behavior of entropy in N- 
body systems with constant total energy Eq, tq is the relevant parameter. It is ro that governs the clumping process 
within a bounded cluster of particles. In particular, ro dictates the extent to which particles can form binaries. This 
is verified by a number of A-body simulations which have been performed for scenarios ranging from formation of 
star clusters to clusters of superclusters |^0|-|l2|. Invariably, the softening length sets the degree of clumping that is 
observed: As the softening length ro is decreased, particles clump tighter and fall deeper into the central high density 
cores. 

One can relate this behavior to the arrow-of-time question as follows: In N-body codes the softening length and 
the mesh discretization scale are equivalent insofar as that, below either, no clumping takes place. (The gravitational 
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force tends to zero and we lose all clumping information.) Hence, by increasing ro, the universe becomes effectively 
more homogeneous, and gravitational entropy decreases. In other words, by changing the mesh size, we change the 
entropy. This is an example of "coarse graining." 

Furthermore, in A^-body simulations, some particles migrate to a dense central core at the expense of ejecting a 
few from the system. Thus, the minimum interparticle separation is a decreasing function of time. We would then 
expect the minimum separation to enter into time-dependent limits of the coordinate integrations. As the separation 
decreases, the depth of the potential well increases and the phase space also. In terms of the two-particle model, if we 
associate ro with the minimum interparticle separation, then as ro decreases, the absolute value between X2 min and 
X2 max in the limits of integration would increase and as well (Fig. |l|) . 

To sum up, the "air-track" model is useful in that it exactly illustrates the calculational procedure we will follow, and 
by reasonable interpretation of r it correctly predicts the behavior of gravitational entropy in A^-body simulations. 



III. HARMONIC OSCILLATOR 



Consider a 1-dimensional system of A" simple harmonic oscillators with Hamiltonian 

JY 



H =2 

»=i 



A N 



and potential 

n&) = §£#- (9) 

With the replacement A~ — ► 3N this may be regarded as a 3-dimensional system with N oscillators in each of the 
x, y, z directions. The phase space for this system can be evaluated analytically. With the canonical momenta 7Tj = fa 
we have, for any form of potential, 

d<pN / dfa / dirN / diri, (10) 

J J -In J-li 

with £^ = 2(H — V) — YliLn+i n i- ^ we further note that £^ = i^+l ~ then each integral is of the form 

' 71 + 1 

i 2 
r( n±2 } 

where x n = n^/i^- (This integral is the so-called Beta function fHH ), 

The final result after A~ integrations over the momentum variables yields the volume of an A^-sphere of radius 
y/H-V and we have for the remaining coordinate integrations 



^ = 7^+27/ Wn-. I dfa{H-V) N '\ (12) 



WW) 

where now l 2 n — 2H/X — 5Zi=«+i 4%- I n ^he case °f the harmonic oscillator, the (^-integrals are identical in form to 
the 7r-integrals. Hence, one merely continues the procedure another N times, and noting that = 2H/X, one gets, 
finally 

A^/2r(7v + i)' [6) 

The same result can be obtained more simply using Af-dimcnsional spherical coordinates. However, the method 
outlined here can be applied to a more general class of potentials. 

Because Eq. (|l3|) will prove central to much of our analysis, it is worth convincing ourselves that the result is 
meaningful. We first note that decreases as A increases, in accord with our notion that a stiffer spring constant 
confines the oscillators to a smaller region of phase space. Note also that A — > =>■ ft — > oo. This behavior is 
equivalent to that of classical free particles in an infinitely sized box. Indeed, by setting V = in (|l2|), performing 
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each (^-integral over the volume of the container and letting 2 N I 2 — > (2m) N / 2 , for massive particles, one can recover 
the usual expression for the entropy of an ideal gas. 

As a further check on Eq. (|l3|), we point out that if one makes the identification H = E = N/ ((3) with j3 = 1/kT, 
then S = k In fl agrees in the classical limit with Einstein's formula for the entropy of N harmonic oscillators S = 
kN(l - ln/3w), where u = V\ @. 

Finally, the usual definition of phase space is the phase space of a shell around Eq: il(AE) — Sl(Eo < E < Eq + AE). 
This is given by the differential of (|l3[). From this differential one can derive the partition function 

Z = jy { AE)e-^E _(*)". ,14) 

Although this formula is not found in texts, if one calculates Z for N oscillators with the ^-function approach textbooks 
apply to free particles jljj, one arrives at the same result. 

With these checks it appears that Eq. ( [i"3| ) gives a reasonable and meaningful expression for the phase space 
of N harmonic oscillators. Perhaps the most important (and useful) feature of (|l3|) is that it merely considers the 
amplitudes of the fa. It ignores the phases. In fact, for its interpretation as phase space we must assume random 
phases for the oscillators. Without this assumption, the motion of the system cannot be considered ergodic and the 
entropy is not defined. Effectively we are regarding the oscillators as a microcanonical ensemble, in which one does 
not know the exact energy distribution. However, one could use the definition S = — ^pilnpi, which to a high 
approximation is equivalent to S = lnfi, and apply it to other distributions as well. 



IV. NEAREST NEIGHBOR POTENTIAL 



The technique used to derive Eq. (|l^) can be used without modification for other potentials of the form V ~ <j) a for 
even powers of a > 0. In addition, as an important application to our analysis, we consider the "nearest-neighbor" 
potential with Hamiltonian 



1 N \ N 



(15) 



i=i 



i=l 



Note that the product — fa+i) is a discrete approximation to the gradient dip/dx; the spatial scale of the 

gradient is set by l/y/X. 

With the substitution r\i = fa — fa+i, the phase space for this nearest neighbor potential can be evaluated in the 
same way as the harmonic oscillator. After integrations over the fa and N-l integrations over rji the result is 



CI 



(2nH) 



N-l/2 



A(JV-i)/2r(AT + l/2) 



(16) 



The lower dimensionality of fl arises from the fact that the last f] used in deriving (|l6| ) is r/jy-i — fay-% — 4>n- The 
final integration, however, requires one to specify boundary conditions to ensure the dimensionality of momentum 
space equals that of coordinate space. A natural choice is periodic boundary conditions, such that fay+i = fa, then 
Vn = 4>n — fa- The assumption of periodic boundary conditions adds an extra (cf>N ~ <fii) 2 to the first integral, which 
can be handled by "completing the square" and pushing the unwanted terms up to successively higher integrals. 
However, a rather tedious calculation shows that, surprisingly, the extra terms vanish after N — 1 integrations. In 
other words, i^-i = 2-ff/A. The dimension has not increased. One therefore is still required to specify limits on fay, 
which we take to be ±y/2H/ A. Then, for periodic boundary conditions, 



N 



f^N \ N ' 2 T(N + 1/2) 



(17) 



However, we note that one might instead impose "free-floating" boundary conditions such that (j)N+i = <Pn, and 
merely specify that the limits on fay are ±^2H/\. In this case the result is the same as ( |l7| ) but without the y/N 
in the denominator. When logarithms are taken, both results are identical to the harmonic oscillator case except for 
insignificant numerical factors. 
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V. APPLICATION TO THE ELECTROMAGNETIC FIELD 



As a sample problem whose technique will carry over to the gravitational case, we now apply these results to the 
electromagnetic (EM) field. Because the EM field can be modeled as a collection of harmonic oscillators, we expect 
the phase space to reflect Eq. ([l3|). To show this is the case we assume a constrained Hamiltonian of the form 16 

H = ±J(£ 2 +B 2 )d 3 x, (18) 

where £ and B are the electric and magnetic field densities. It is important to remember that in the Hamiltonian 
formalism, the canonical variables are not the densities but the full field quantities; in this case ir = E and q = A, 
the vector potential. 

To write Eq. (|18|) in terms of E and A, we first discretize H as follows: 



N N N 



2^^^ N x N y N z 

x y z * 



Here, the sums are understood to be over the x, y and z coordinates covering the enclosed volume L x L y L z . We have 
also approximated dx as LAN/N, with N being the number of oscillators in each direction and AiV = 1. Let us 
further restrict our attention to transverse waves propagating in the z-direction. Then the x and y summations can 
be easily evaluated with the result 

1 N r 3 
H= 1 -N 2 J2(£ 2 +>3 2 )^, (20) 



N 3 

where L 3 = L x L y L z and N 3 = N x N y N z . 

Now, £ ~ y / e/L 3 , if e = energy. Similarly, if A is the potential density, then B ~ ^Je/L 3 ~Vxi~ A/L. Hence 
£A~e/L 2 . 

The product of the canonical variables, EA must equal an action = eL in these units. Consequently, 

£AL 3 = action = EA, (21) 
and the proper scaling becomes E = £{L/N) 3 / 2 and B — B{L/N) 3 / 2 . The Hamiltonian ( pfj| ) can then be written as 

JV 

2 



H = W{E 2 + B 2 ), (22) 



where H = H/N 2 . 

To evaluate the phase space below this Hamiltonian, we note that for transverse waves, V = V z ; E z = B z = 0; 
and B = n x E, where n is the propagation vector. Faraday's law, V x (E + A) = 0, implies that A z — and 
B = —iA y z + jA x z . Thus the Hamiltonian can be approximated as 



^ 1 N 



H=-J2(, E U*)+E 2 y W)+V(A), (23) 



2 

i—i 



with potential 



A N 



V ^ = 2 E [(4,(0 ~ Mi + !)) 2 + (4.(0 -Mi+ I)) 2 ] ■ (24) 
i=\ 

and where A sets the spatial scale. We see that in this approximation, V is just given by the nearest-neighbor 
potential, calculated in §IV, with the phase space given by Eq. ([l7]). In this problem, however, the phase space is AN 
dimensions, hence substituting 2N for N in Eq. ( [L7] ) yields 

HT~ (2n~H) 2N 

^ EM = i^Wr(2NTWy (25) 



G 



For the interpretation of fl as phase space we need to assume the phases of the electromagnetic waves are random, 
which merely means the source is incoherent. This is, in fact, the general case. 

A closed solution to the problem requires an evaluation of N, the number of oscillators. Because we are primarily 
interested in the time dependence of f2 (which is here time independent), it is enough to know that N is finite; in the 
quantum limit it will be the number of photons. 

To make contact, however, with the usual expression for the entropy of electromagnetic radiation, we imagine 
transverse waves in a three- dimensional box, assuming each direction is independent. The phase space for that system 
is obtained by letting TV — > 3N in the above equation. We assume that H = 3Ne , where e = lu/(2tt) is the average 
ener gy p er oscillator and the coupling constant A = uj 2 . With Stirling's approximation T(N + 1/2) as y/2TrNe~ N N N , 
Eq. (|25| ) yields S = In 51 ss 6iV(l — In 2). For a photon gas at temperature T, the energy of most photons is of order 
u> = k ~ T, where k is the wave vector. In three dimensions, the number of states is proportional to the volume in k 
space for a sphere of radius The mean number of photons at a temperature T is thus proportional to N ~ k 3 ~ T 3 , 
and we recover the usual scaling for the entropy S ~ N ~ T 3 . 

We also point out that ( |2~i| ) gives some insight into the question of coarse graining. The concept of entropy is 
subjective in the sense that to calculate an entropy requires that an averaging procedure be selected. If one regards 
the coupling constant in ( pi] ) to be A = N 2 /L 2 , where L is an arbitrary length scale, then by increasing L, one 
increases the volume per oscillator and hence increases the phase space, as can be seen from (|2^). Therefore the 
coarse graining scale evidently appears in these calculations as the coupling constant. 



VI. EXTENSION TO GRAVITATIONAL WAVES 



The extension of the previous formalism to gravitational waves is fairly straightforward except for one crucial point, 
which we discuss below. 

We first consider inhomogeneous perturbations of the spatially flat metric 



ds 2 = a 2 (rj) —drj 2 + (Sij + hij(r], z)jdx l dx^ 



(26) 



where rj is the conformal time, a(rf) is the expansion scale factor and the hij <C 5ij represent gravitational wave 
perturbations. Their equation of motion can be found by expanding the Einstein action to second order in the 
perturbation variables h^ u . The result is 



1 

64^ 



a 2 (h 2 -h u )d 4 x, 



(27) 



where (•) = djdr\ and (') = d/dz. Eq. ( |27| ) is the action appropriate for singly polarized gravitational waves in 
the transverse traceless gauge. The variable h (= h xx = —h vy ) represents the single degree of freedom for the + 
polarization state. (Brandenberger et al. || have shown that a similar form is achieved even when one considers two 
polarizations.) 

We will find it convenient (particularly when making a connection to density perturbations) to work with a trans- 
formed perturbation function <fr = a/i/v327r. The Lagrangian density can then be written as 



a2 



By definition, if <p is the canonical coordinate, then 7r = dC/dcj) = 
found to be 



1 

H =2 



A'2 



(28) 

and the Hamiltonian density TL = irq — C is 

(29) 



Following the same procedure used in the EM case, we find for the Hamiltonian 

N 



a— 2\ 
-<Pi i ' 



(30) 



where W = n(L/N) 3 / 2 , = 0(L/7V) 3 / 2 , H = H/N 2 , and H = fHd 3 x. We see that the Hamiltonian contains 
potentials similar to the others we have considered with one crucial difference: the sign on the last term. In the 
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matter dominated period, a ~ r/ 2 and a/a > 0, hence the sign on the quadratic term in Eq. ( |30| ) is negative and we 
have a nearest-neighbor potential plus an "antiharmonic oscillator" (or inverted) potential. 

The inverted nature of this term is due to the background curvature of spacetime and its rate of expansion. The 
potential, then, serves as a reflection barrier in an unbounded phase space. Any calculation must therefore include an 



arbitrary cutoff. We discuss this point in detail in section § VIII. There we show that Q can be calculated by the use 



of hypcrgcomctric functions with a result that is formally similar to that already achieved for the harmonic oscillator 
and we can continue to use Q ~ H N /\ N / 2 to compute the time dependence of f2. 

To compute £l(rj) we imagine that fl is constant on each hypersurface of constant time. Thus a/a can be taken as 
A, the coupling constant. We then need to compute H(rj) and X(rf). To find H{rj) note that the equation of motion 
for h resulting from varying the action (E7f) is 



h + 2-h-h" = 0. (31) 
a 

Because the waves are linear perturbations, they do not interact except through a linear superposition. The time 
development of each individual component (or a wave of a particular frequency or phase) evolves according to Eqn. 
(|3T|), which can be regarded as representing a family of solutions. That is, we may assume a separable solution 

^2 Chj + —kj -h'A = (h+—h- h") &iai = °> ( 32 ) 

j 3 

with arbitrary or random phases ay, where here j is an index over the different waves (not over z). This is, in effect, 
saying that different spatial regions are taken to be oscillating independently of one another, or that the source is 
incoherent. We therefore assume the perturbations to be random, and that the field variables describe, not a singly 
polarized wave, but an ensemble of incoherent plane waves. Entropy is thus attributed to the lack of knowledge in 
the exact field configuration. 

With a = a a rf for the matter dominated period and assuming h ~ e zkz , Eq. ( |3l|) has the solutions 

h<xrr Z/i J±*/iW)e ikz , (33) 

where J±3/2 are Bessel functions. To construct the Hamiltonian (^p|), we then sum over the coordinate z. Simplified 
expressions can be obtained by substituting the standard asymptotic {krj <C 1 and krj 3> 1) forms of J±3/2- We can 
then write for the metric perturbations in the limit krj <C 1 

h= \h l {kr ) )- 3 + h 2 \e ikz , (34) 

where hi and h 2 are constants and </> = ah oc (ft-if? 1 + h 2 rj 2 )e tkz . The constants hi and h 2 can thus be interpreted 
as defining the decaying and growing mode solutions respectively. In this limit spatial gradients are negligible. For 
krj 3> 1, we have 



h oc y ( — j cos(ki]) +sin(/c7y) e lkz oc (kr])~ 2 x [oscillations]. (35) 

and 4> oc constant x [oscillations]. We note that fc)j > 1 {krj <C 1) represents perturbations with wavelengths much 
shorter (longer) than the Hubble radius (usually referred to as the "horizon" ) . 

The Hamiltonian (|30|), in the limit krj 3> 1, then becomes H oc n 2 + k 2 <f> 2 which is simply the harmonic oscillator 
Hamiltonian at a fixed time with coupling constant k. H therefore oscillates in time at constant amplitude and we 
have for the phase space 

Ct oc ^ N / 2 x constant x [oscillations] . (36) 

As expected in this approximation, the phase space does not change. Recall that H is defined on a single time slice. 
However, assuming incoherency in time as well as in space, one can average H over several cycles by defining a general 
4-Hamiltonian 

W« = J d 4 x[n 2 +k 2 4> 2 -^]. (37) 
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Over intervals of time greater than the dynamical time, this will be a monotonic function and, in the krj 3> 1 case, f2 
will be strictly constant. However, in a nonlinear regime, Q would increase, which is encouraging for the interpretation 
of In f2 as entropy. 

For k-q -c 1 spatial gradients are, again, negligible and we have H oc ir 2 — a <j> 2 /a and oc H N /{d/a) N l 2 . Therefore 

H oc ( ^ 4 ' and Q oc ( ^ ' !° r f "^ m ° d f ' (38) 
I r} , I ?7 , for decaying modes. 



with the caveat that we have not yet shown (see § |VII]| ) that for the inverted oscillator this form of Vt is justified. 

At first these last results strike one as strange because as seen from (pij), for growing modes h is frozen- in at 
superhorizon scales. That is, there are no oscillations and the assumption of random phases is not well motivated. 
In that case, the phase space trajectories are known precisely and f2 = 0. One can see this clearly by examining the 
Hamiltonian in the variable h. For superhorizon growing modes, this Hamiltonian is zero, and the increase in O above 
is entirely due to the expansion of the universe (i.e., 4> = ah + ah — ah). Only the a term causes <j> to increase in 
amplitude, and hence increases the effective coarse graining scale and f2, in accord with the gravitational arrow of 
time. (The decaying modes on superhorizon scales also no longer oscillate but damp out monotonically at a rate faster 
than the universe expands; the asso ciated entropy thus decreases.) We will discuss the significance of the growing and 



decaying modes further in § VIII B . For now we point out that, certainly for the growing modes, Q is nonzero only if 
one continues to regard the phases as random. Otherwise, if the phases are assumed known, then the entropy is zero 
(or constant), in agreement with Brandenberger et al. |8j. It would be of interest to establish a more quantitative 
comparison of our results to Brandenberger et al. 

These considerations suggest that our definition of entropy is only appropriate for subhorizon scales. This may 
actually be an advantage, because in order to ensure that fl is finite, we must ensure that the number of modes N 
must also be finite. This requires us to put the system in a box and consider only a finite spatial region. The horizon 
thus provides a natural upper limit to wavelengths. An absolute lower limit can obviously be chosen as the Planck 
scale. We will find that similar considerations are necessary for radiation perturbations (below). 

VII. DENSITY PERTURBATIONS 

The analysis of the previous section can be repeated for density perturbations in dust- and radiation-filled models. 
In the longitudinal gauge, the spacetime metric is written as 

ds 2 = a 2 (j]) [-(1 + 2$(i7, z))dj] 2 + (1 - 2$(t7, z)) lrj dx 1 dx J 1 , (39) 

where 

K 



1 + j (x 2 + y 2 + z 2 ) 



(40) 



<f> is the gauge invariant gravitational potential, and K, = 0, -1, +1 for flat, open and closed universes respectively. 
Mukhanov et al. [[l7) give the following general equation for adiabatic density perturbations: 

U - c 2 s u" - -u = 0, (41) 

where 

' i<!> f 2 a 4-~) , (42) 



V47T V a 



3 a I a a x 



2a 



2-5 - - , (43) 



and (•) = d/drj, (') = d/dz and c 2 is 1/3 for radiation and zero for dust. The corresponding action from which the 
(ADM) Hamiltonian and equations of motion are derived is given by writing the Einstein action, 

I = / R^-gd^x - ( eV=Q~d 4 x, (44) 



16ttG . 

where e is the energy density of matter, in terms of the ADM metric and expanding to second perturbative order 
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A. K, = 0, flat universe 



1. <? a — 1/3, radiation 



For radiation, a ~ 77 and 0/0 = 2/?? 2 . Assuming the spatial form ~ e , the field equation pH ) becomes 

2 A; 2 



— u+— u = 0. (45) 



The general solution to ( f45| ) involves Bessel functions similar to Eq (|3^). The asymptotic super horizon [krj « 1) 
solutions are 

u=(u 1 ri 2 + u 2 ri- 1 )e ikz ) (46) 

representing both growing and decaying modes. As in the case for gravitational waves, these solutions are taken to be 
a family of functions with random phase angles a. From now on the presence of these phase angles is understood but 
we do not write them out explicitly. In addition, we note that the results presented here are independent of the exact 
form of perturbations, and we could replace e^ fcz ) Ylj e%aj with an arbitrary function of the three spatial coordinates. 
For krj^> I 

u cx (cos -p: + sin ) e ikz , (47) 

and, as for gravitational waves, perturbations on these subhorizon scales are oscillatory. 

The Hamiltonian for the case kr\ -C 1 is H cx 7r 2 — 8u 2 /8 which results in the following evolution 

rj f V 2 : j o \ rf N , for growing modes, ,. Q , 

H cx { _a and il cx { _ 3N c ° . , (48) 

\ T) , [ 77 , for decaying modes. 

For kr\ ^> 1, we have H cx 7r 2 +c 2 w' 2 which oscillates at constant amplitude, and therefore is constant over sufficiently 
long intervals of time. Notice that these results for radiation perturbations are identical to those of gravitational waves 
and the remarks concerning the superhorizon application of our definition of entropy apply here as well. 



2. d 



0, dust 



In this case a = a^rf and 6/9 = Q/if. Then (|l|) becomes 



with solution 



or equivalently 



6 

-^u = 0, 
T 



(u lV 3 + U2TJ- 2 ) e l(kz \ 



$ = (hi +u 2 rj 5 ) e 



(kz) 



(49) 



(50) 



(51) 



where u\, u 2 , ui and u 2 are constants. 

Notice the important point that Eqn. (^) does not suggest a natural scale (the horizon in particular) for modes 
to grow or decay as found in the gravitational wave and radiation cases. However, the horizon scale does appear in 
the expression for the density fluctuations Sp/p (l?]] 



6_p 
P 



-(fcV + 12)ui - (fcV - 18>r 5 u 2 



(52) 



Thus, in distinction to the previous cases, one should not examine <f> to determine whether the modes are frozen-in or 
not. One should rather examine (|52|). We see that on subhorizon scales (ki] 1) Sp/p exhibits growing modes, i.e., 
actually collapse takes place even while $ remains constant. On scales larger than the horizon, 6p/p remains constant 
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for the dominant growing modes. Given that dp/ p is constant on superhorizon scales for growing modes, this once 
again suggests that our definition of entropy should be restricted to subhorizon regimes, consistent with our earlier 
results. 

For K, = dust we have simply H oc n 2 — 8u 2 /8. From ( |50|) and the results from § VIII B for the inverted oscillator 
potential, we find 

H oc ( t 6 ' and Q oc ( ^ ' !° r f ° wing mod f' (53) 
[ rj , I >7 > i° r decaying modes. 

We again point out that the fact that H and £1 are time-dependent while the conformal metric components (or 
equivalently $) are constant is not a contradiction. The growth of H and Vl is an indication that collapse is taking 
place on some (subhorizon) scale. Note once more the monotonic behavior of these quantities. 

B. K, = ± 1, dust-filled open and closed universes 

The gauge invariant potential in the case of a dust-filled open universe can be written as |l7j ] 

2 sinh 2 77 — 677 sinh 77 + 8 cosh r/ — 8 sinh 77 

$ = (cosh, -1)3 + ^ (cosh, -1)3 ^ 



where a ~ (cosh 77 — 1) is the expansion factor, c\ = uie^ kz+a ^ and C2 = U2e ltykz+a \ Also 

(cosh??- 1)V 2 $ + 9$ - 6cij . (55) 



5p _ 1 

7~ 3 



Expanding Eqn. ( |54| ) in the small time limit 77 <C 1, we obtain the flat space solution (pl|). Eqns. ( |43j ) and (p2[), taken 
together with the approximate asymptotic solution for the scale factor a ~ t? 2 , yields the same result for H and f2 
as the flat space case (p3|). In the opposite late time limit, 77 ^> 1, the hyperbolic functions become exponentials and 
$ oc cie~ n + c-2e~ 2ri . As expected, $ decays in this limit, and ([55]), along with the asymptotic form of the scale factor 
a ~ e n , shows that the matter density fluctuations do not grow on either super- or sub-horizon scales. Eqs. ( f43| ) and 
([l2|) yield 9/6 = constant and u oc ci + C2e~ T> . H and 17 then evolve as 

( constant , , / cons ^ an ^ 1 ^ or "growing" modes, . 

1 e _2,) , ' 1 e~ 2Nr i , for decaying modes. ^ ' 

This is consistent with the fact that Sp/ p — constant for the dominant modes and no collapse takes place on either 
sub- or super-horizon scales. 

The growing modes for the closed model (K. = +1) can be obtained by letting 77 — * 777 in (p4|). In this case 
< 77 < 2-7T. For 77 <C 1, the closed model gives the same result as the open and flat cases. Eq. ( j55[ ) also holds for 
the closed model, and so the same consistency among the behaviors of H , and dp/ p is found here as well. There is 
no asymptotic limit 77 3> 1 in the K = +1 case. However, by expanding ( |54| ) around 77 = n, we find to lowest order 
that $ oc 77 ~ constant; 6/8 ~ constant and that therefore H and f2 are constant as well, again consistent with the 
density perturbations Sp/ p ~ constant. In the neighborhood of maximum expansion, then, the model acts like the 
open case. As recollapse takes place one finds that toward the singularity 77 ~ 2ir, $ oc e -5 , u oc e -2 and 8/6 oc e~ 2 , 
where e = 2ir — rj. These dependences yield H oc e~ e and oc e~ 5N . Since e is decreasing, all these quantities are 
increasing, as expected; 5p/p grows again as well. Taken with the behavior at 77 ~ and 77 ~ it we see that H and CI 
behave monotonically as required. 

VIII. ANTIHARMONIC OSCILLATOR POTENTIAL 
A. Time Dependence of Q, 

We now consider in detail a Hamiltonian of the form found in Eq (|3^) without the gradient term, that is, 

N N 

J ff = A^x 2 -^y 2 . (57) 

i=l i=l 
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We first justify the exp ression s found for f2 in §VI, then interpret the meanings of H < and H > for the inverted 
oscillator potential in § VIII B . The inverted nature of the potential in (^7|) results in a reflection barrier and a phase 
space that is unbounded. To compute f2, then, we will need to put in arbitrary cutoffs to the allowable phase space. 
How this is done will become clearer below. 

First consider the case when H > 0. From (|57]) the yi then correspond to the canonical coordinates <pi and the Xi 
correspond to the previous 7Tj. One can easily evaluate O using ^-dimensional spherical coordinates. Integration over 
the yi coordinates yields 



and then over the 



n = 



Q,y — 



TT N N 



r N/2 



r( 



N±2t 
2 I 



N/2 



' N 

.i=i 



H 

T 



N/2 



r( jV±2)r(JV±2) 







' A 2 " 


(i: 








' JR 2 >H/\ 


H R - 1 



-i N/2 



R N ~ 1 dR, 



(58) 



(59) 



where R 2 = J2iLi x i- 

This integral will be unbounded as R 2 



We therefore let R 2 



maximum(^]^ 1 x 2 ) be the assumed cutoff in 



x-space (which here corresponds to momentum space). Further defining u Tl 
the above expression becomes 



R max^/ H and w = 



N ttN 



n 



TT N H 



I) 



N/2 



(it, 



r(f )r(*±*) V^a; 

where J 7 is a hypergeometric function 

'2-N 2 + N 4 + TV 



1) 



(JV+2)/2 



J" 



2-iV 2 + iV 4 + iV 



[l + (t 



•.s i(iV-2)/2 TV/2 7 

l)iuj v " u> ' aw. 



2 ' 2 ' 2 

The function T is absolutely convergent for \u max — 1| < 1. In that limit, the power series |l3j for T gives 

2n N H (N-2)/2 A 



n 



(7v + 2) r(f)r(^)e JV/ 



rJ2 



JV+2 



2 max i 



(60) 



(61) 



(62) 



which can also be obtained by direct integration of (|6C| ) if one keeps only the tu^/ 2 term in the integrand. Thus to 
evaluate the time dependence of Q we will consider 



n 



-1 I'.-vyi , , ■>■ 



az jt(N-2)/2 



for 



^Caz = constant 



k N/2 

where k = 2£ is the "spring constant" . With the definition of u maa; this can be rewritten as 

R N 



k N/2 



for 



1 



1- 



(63) 



(64) 



assuming u max is constant. That is, which form used depends on which variable is assumed constant. Conceptually, 
it is easier to visualize the meaning of R m ax, which puts an absolute limit on the allowable momentum of oscillations. 
We stress that these limits are meant to be constant in time. If, however, we allow Rmax to evolve with H , then we 
get the condition u rnax = constant in time. The parameter u max effectively scales the ratio of the allowable kinetic 
energy to the total energy H . 

We can also derive an analogous scaling in the opposite limit, R max , u m ax ~ * °°- hi this case, Eq. (p0[ ) is 
approximately 



^N 



n 



2T{ N±2 )r{ N±2 ) ^ 



N/2 R 2N H N ii N 

p2N ^max _ a max 

U max CX- k N/2 °^ foN/2 ' 



(65) 



Now we turn to H < 0. In this case the meanings of the Xi and yi in ( |57| ) are reversed. If we let H = \H\, then a 
repetition of the previous analysis gives 
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^max 



/k N/2 , u max = constant, 



~ { kHl»-WR£g, R 2 max = constant ~ H/k, (66) 
k N/2 R™„, i£L„ = constant — > oo, 



max 



where the "spring constant" is now fc = 2A and R ma x corresponds to a cutoff in </>-space (or places a limit on the 
amplitude of oscillations). 

At this point we reiterate that our approach is to compute phase space with the assumed cutoffs at each time 
slice and then allow the system to evolve in time. Since the Hamiltonian is, in general, a function of time, it is 
reasonable to impose cutoffs that scale with H. This implies that we should hold u max constant, thereby preserving 
the self-similarity in the energy distribution. This choice of cutoff is also computationally convenient in that it results 
in a scaling for £1 that is similar to the harmonic oscillator case, namely fl oc H N /k N / 2 . However, we have verified 
that the other cutoff criteria (constant i? 2 rlQ a:) gives qualitatively the same behavior as the constant u rnax case for 
both growing and decaying modes. 



B. Meaning of H > and H < 

In the previous section we considered fl for both H < and H > 0. We now wish to explore the meaning of 
positive and negative Hamiltonians in the present context. Classically, one associates negative energy states with 
bound systems and positive energy states with unbound systems. Here, however, the situation is slightly different. 

The equations describing the evolution of dust and supcrhorizon radiation and gravitational wave (with the iden- 
tification u = <f> = ah) perturbations can be unified into a single differential equation u — cur\~ 2 = 0, where c = 6 
for dust and 2 for radiation and gravitational waves. For both cases we can write the solution as u = Arf 11 + Br/ 112 , 
where A and B are constants and (n\, n^) — (3, -2) for dust and (2, -1) for radiation and waves. In short, A and 
B define the growing and decaying modes respectively. From our previous solutions to the equations of motion on 
superhorizon scales, the Hamiltonian H ~ u 2 — cu 2 r]~ 2 can be written as 

_ f 3AV - 2B 2 ?r 6 - 24AB7T 1 , for dust , . 

y 2A 2 i] 2 — i? 2 ?7~ 4 — SABrj^ 1 , for radiation & waves ^ ' 

Now, note that A = results in 

H _ f -2B 2 J]- 6 < 0, for dust 

1 —B 2 r]~ 4 < 0, for radiation & waves ^ ' 

That is, in both cases, H < corresponds to a decaying mode. Similarly, setting B = leads to 

H _ f 3AV > 0, for dust 

1 2A 2 r) 2 > 0, for radiation & waves ^ ' 

In other words, H > corresponds to a growing mode. 

It is important now to attach a physical picture to these results because they are, in a sense, opposite from what one 
intuitively expects from a particle model. In a particle model, one associates H < with bound systems undergoing 
gravitational collapse. Growing modes, then, correspond to H < and particles moving together. 

However, it is crucial to bear in mind that we are considering not a particle model but an oscillator model, where 
growing modes correspond to increasing amplitudes of oscillation. One therefore can imagine a lattice of points 
undergoing perturbations that eventually lead to gravitational collapse. As the perturbations grow, the grid points 
move further from their initial unperturbed, or "uniformly" arranged positions. For decaying modes, the grid points 



relax to their homogeneously spaced positions. This is why in §VI, Q grew for growing modes and decreased for 
decaying modes. In the oscillator picture, then, increasing inhomogeneity automatically gives an increase in phase 
space and hence gravitational entropy. 

We can also make contact with the "qualitative cosmology" approach of Hamiltonian cosmology. The turning 
points of trajectories with the inverted potential Hamiltonian take place when the momenta are zero and H = V. 
Since for the inverted potential V < 0, H > necessarily implies H > V. The motion here is "unbounded," in the 
sense that there are no turning points and perturbations continue to grow. For H < 0, we have \H\ — CX^i^i 
at the turning points. The potential barrier is thus an N-sphere of radius r = y/\H\/£. However, in general, 

J2iLi 4>1 — \H\/€ + ^Ei=i "f/C > \H\/€> so tne wor ld point is actually outside this sphere. 

For decaying modes, the sphere shrinks in time and the world point attempts to catch up with it. However by 
comparing u for decaying dust and radiation modes with the time dependence of r for the potential barrier, one 
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easily shows that the system point can never catch up with the barrier in finite time. The barrier, then, serves as 
an attractor for the decaying modes but it is never actually reached except in an asymptotic sense. This picture is 
similar to that of the Bianchi cosmologies, in which the universe is often represented as a point moving in a potential 
well. We turn to Bianchi IX cosmologies now. 



IX. BIANCHI IX COSMOLOGY 



For the Bianchi type IX cosmological models, we adopt a metric of the form 

ds 2 = -dt 2 + e 2a (e^)..aV, (70) 

where a — e a is the mean expansion scale factor, a 1 are the dual 1-forms for the rotation group SO(3,R), and 
(e is an exponential of a 3 x 3 symmetric traceless matrix defining the anisotropy of the spatial hypersurfaces 
and parameterized as 

= diag + V3/3_, (3+ - V3/3_, - 2/3+||. (71) 

The Bianchi models are anisotropic but homogeneous cosmologies, so by definition they cannot show the effects of 
gravitational clumping. Nevertheless, there are three reasons for investigating Type IX. First, it can be conveniently 
cast into a Hamiltonian form and a phase space can be formally calculated. Second, if one regards anisotropy as the 
long- wavelength limit of inhomogeneity, we might hope make contact with our previous results. Finally, it provides a 
transition to the full ADM formalism, which one will necessarily employ in nonperturbative models. 
The ADM Hamiltonian for Bianchi IX is 

H 2 = p\ +p 2 _ + 36Tr 2 e 4a {V (J3+, /3_) - 1) . (72) 

Hence, Bianchi IX can be cast into a system of two degrees of freedom (meaning two canonical pairs.) In jT^ ) 
V is Misner's anisotropy potential, which is a function of the canonical coordinates /3+ and /3_, the independent 
components of the metric anisotropy. The precise form of V is: 

V = 1 + \e- &f3 + - \e~ 2f3 + cosh2V3/3_ + |e 4/3+ (cosh4 v / 3/3_ - 1). (73) 

This potential (shown in Fig. |^) is symmetric about the /3_ axis and has exponentially steep walls. For large 
isocontours of V (> 1), the potential exhibits a strong triangular symmetry with three narrow channels that extend 
to infinity. For V < 1, the potential is closed and asymptotic (/3± <C 1) isocontours describe a circle. The motion of 
the universe point in this potential well is chaotic p0| , so we can regard any region of phase space to be filled with 
equal probability and the concept of an associated entropy is reasonable. 

The phase space is formally calculated as we have done with the other cases: 

Mix = f d P+ J d P- J d P+ J d P- ( 74 ) 

To facilitate integration, however, we eliminate p- in favor of H. The integral then becomes from ([F: 

fH max ,/?+,_ f (3- ma „ f +e d 

n IX = J HdH / d/3+ / df3- / , (< r >! 



0+min J - 1 \jV 2 -y\ 

where £ 2 = H 2 — 367r 2 e 4a (^ — 1). We note that this problem is very similar to the two-particle model of §0, except for 
the more complicated potential. The integral over p+ is simply an arcsin and the result after applying the boundary 
conditions is tt. The lower limit on H is 0, and due to the symmetry of the potential, we can take the /3_ limits to be 
and 

—max i the maximum value of /3_, and double the result. Therefore we are left with 

fflmai fP-max Pp+max 

n IX = 2Ti / HdH / d/3_ / d/3+. (76) 

JO JO JP+min 

The remaining integrals are evaluated numerically. To do this requires first determining the limits of integration. 



As with the two-particle model, we set limits by equating the momenta in (72) to zero and finding the reflection 
points. In other words, we demand that H 2 always remain positive: 
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H 2 > m^e^iV -l)^V< + 1. (77) 

For a fixed value of H and /3_, we march across the potential well varying (3 + until this inequality is violated. We 
then perform interpolations at the two endpoints to find the minimum and maximum values of /3+. Then /3_ is 
incremented and the process is repeated. The limits on /3_ are found in a similar manner. For large values of /3±, we 
treat the equipotentials as equilateral triangles. For intermediate values of (3± we take into account the deformation 
of the contours and follow them part way into the channels (the area here becomes vanishingly small). For a closed 

potential in which [3± <C 1, the area is approximated as a circle of radius \J 0+ + 0L- At each time step the integrals 

in equation ( |76| ) are evaluated with a 40-point Gaussian quadrature scheme. 

To evolve the system in time, we integrate the evolution equations for f3± and a 

h = -W+-\<T**§- + , (78) 

/L = -3d/3_ - \^ 2a ^~, (79) 

Pl+pl~\e- 2a (\-V)] 1/ \ (80) 

using a 4th-order Runge-Kutta scheme. The Hamiltonian is then updated at each time step by 

H = 12-Kae 3a . (81) 



This value of H is then used for the upper limit of the outer integral in ( |76[ ) . 

We note that the /3-integrals basically give the area of the triangular potential. Thus we can estimate the size of 
the phase space as 

n IX = 2-k J H'dH 1 J A w 2tp^, (82) 

where A is the area of largest triangle and the factor of 1/3 is introduced to approximate the size of an average triangle 
in the inverted pyramid of the potential well. Estimates performed this way typically agree with the computed results 
to within a factor of two or better. 

Results of the numerical integrations are shown in Fig. |^ where we plot the Hamiltonian and the volume of phase 
space as a function of a. We note that the limit a — > — oo corresponds to the "Big Bang" singularity. One of the most 
striking features is that both the Hamiltonian H and the phase space f2 are seen to oscillate. We now demonstrate 
that these oscillations are real. From (|82|) we have Vl ~ AH 2 . Then 

dtt , dH 9 dA 

2AH— + H 2 — , 83 

da da da 

where from (jT^) we have the fundamental equation 

dH 727r V(V-l). (84) 



da H 

Assuming the boundary triangles for the unbounded (open potential) phase space are equilateral, the enclosed area 
is approximately A ~ V3/3 2 . Also, for large values of /?_, the asymptotic form of V is from ( ff3"l) V ~ (l/3)e 4v ^ /3_ . 
Thus A - (ln3V0 2 . Now, at the potential wall, @ shows that V » H 2 e~ 4a /36ir 2 and we can write 



.4 



, , H 2 e~ 4a 
In , 

12tt 2 



dA / l dH \ , [H 2 e- ia \ 



Analytic approximations for the behavior of Q can be found for two limiting cases: i.) "Free-particle trajectories." 
Such trajectories correspond to the plateaus in Fig. ||. In these regions, the universe point is sufficiently far from the 
potential walls that the potential terms in (|72| ) can be neglected. The universe point propagates like a free particle 
with constant "energy" H, so that dH/da ~ 0. 

ii.) "Wall collisions." At or near the potential barriers the momenta in H are negligible so that dH/da ~ 2H or 
H e 2a . In Fig. ||, wall collisions correspond to the places where H suddenly increases. Here the system is gaining 
energy from the gravitational field. During wall collisions the area remains approximately constant so dA/da ~ 0. 
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From equation (p3[) we then have for the free and bounce cases respectively, 



^~-8iJ 2 (21ntf-4a-lnl2^ 2 ), ^ ~ 0, (86) 
da da 

and 

da da 

Equation ( pq ) shows that for large negative a, dVl/da < 0, as observed in Fig. ||. Furthermore this scales roughly 
as ~ H 2 , so as H increases dft/da becomes more negative and the phase space evolves more rapidly, which is also 
observed. 

Note that dfl/da in (^) is positive definite, so at wall collisions always increases and, for a large enough energy, 
increases more rapidly than H. Furthermore, the greater the value of H , the greater the slope, as observed. Together, 
the competing behaviors in the two limiting cases account for the oscillations observed in the figure. 

One should not be disturbed to find oscillations in entropy in this situation. H is time dependent and there is 
no law that says for a time-dependent Hamiltonian, the entropy should monotonically increase. Indeed, the Bianchi 
IX model we have been considering resembles more closely an open system, although the "particles" are not in any 
sense in a canonical distribution, so we cannot compare the size of the entropy fluctuations with those expected for a 
system in contact with a heat bath. 

In terms of findinga monotonic function to call entropy, we also reiterate that this system is homogeneous. However 
two aspects of Fig. || are highly encouraging. For late times (a > 0), the phase space — and hence entropy — is seen 
to increase monotonically in the direction of increasing anisotropy. Furthermore, in the oscillatory regime, along the 
plateaus, where the model most closely resembles a typical "closed" system (E = constant), the entropy is again 
seen to increase in the direction of increased anisotropy. Both behaviors correspond with the notion that anisotropy 
represents the long- wavelength limit of inhomogeneity. 



In the limit of small anisotropy V « 8(/tf + j3 2 _) < 1, Eqns. (78) — (BOj) become 



h + 3-/3± + \p± = 0, 



a 



6? = Pl + fll- e —, (89) 

where we have defined a = e a and a = a/a. Equation (^8|) is similar to ([H]) for gravitational wave perturbations. 
However, the type IX solution is further complicated by Eqn. (|8^) which couples the anisotropy to the expansion 
factor. Nevertheless, at late times, we expect f2 to behave similarly in both the Bianchi IX and gravitational wave 
cases: increase monotonically with increasing anisotropy or inhomogeneity. 



X. COMPARISON WITH C 2 

Penrose had suggested that the square of the Weyl tensor C 2 = C a p S C J " 13 might act as an arrow of time, 
increasing monotonically in time as the universe becomes more inhomogeneous. Of course, this presupposes an initial 
low entropy state at the singularity, in which the matter distribution is homogeneous and the Weyl tensor tends to 
zero. However, Wainwright & Anderson ^ (see also Goode & Wainwright ||) have shown that cosmological models 
which admit an isotropic singularity, contradict Penrose's hypothesis. They also noted that the Ricci tensor diverges, 
but in such a way as to dominate the Weyl tensor. This lead them to propose a weakened form of Penrose's hypothesis 
in which the quantity 

2 ^ ~/S r< a/3 

R 2 ~ R a f)R a P [ ' 

might be the appropriate indicator. However, subsequent work by Bonnor |^] has thrown even this weakened form 
into question. 

Here we calculate the two variants of Penrose's proposal for cosmological density perturbations in an expanding 
flat universe. Assuming, for simplicity, the perturbations to be functions only of conformal time rj and a single spatial 
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coordinate z, the spacetime metric is given by (^) with K = 0. We find (using MathTensor and Mathematica) to 
lowest order in the smallness parameter 1 

r yS r a/3 _ 16($ zz ) 2 , , 

° Q /3 ° 7 <5 - 3^4 ' 

and the solution for <E> is given b y (plf ). During the matter dominated regime, the scale factor evolves as a ~ rj 2 , 
and we immediately see that Eq. (|9 if) does not produce the right behavior for the growing modes. The Weyl tensor 
decreases with increasing time and inhomogeneity. Because the overall time dependence is monotonic, one might think 
to correct this by introducing a negative sign: However, then for the decaying modes C 2 is increasing for decreasing 
inhomogeneity. 

Following the suggestion of Wainwright & Anderson [|| we also calculate 

CjC^ 4a 4 (^ z ) 2 = ± (*.»)' fo 2l 

R af }R af} 9(d 4 - aa 2 h + a 2 d 2 ) 27 ' 1 ' 

Equation ([)2j) does have the correct behavior. In fact, it is interesting to note that to this order the time dependence 
is identical to that found for the corresponding Hamiltonian (|53|), ie. rf and r]~ 6 for the growing and decaying modes 
respectively. 

The generalization to nonflat spacetimes is rather complicated and not qualitatively different from the flat case, so 
we do not include it here. However, we do compute the Weyl tensor for spacetimes of the form ( |26| ) , containing singly 
polarized (+) small amplitude gravitational waves propagating in an expanding universe. In this case 



C^C^ q 4 ( h % z - 4h% z + 2h zz h nv + h 2 J 
R af3 R a f 3 12(a 4 - aa 2 a + a 2 a 2 ) 



(93) 

Noting that R a /3R al3 , to zero perturbative order, is the same as for the metric ([39]), we again find that C 2 alone does 
not produce the right monotonic behavior, but C 2 / R 2 does. To evaluate the latter, we assume the expansion of the 
universe is governed by density perturbations and that the scale factor evolves as a ~ rj 2 . For superhorizon scales, 
krj <C 1, we may ignore spatial gradients so that only the last term in Eq. (|93| ) survives. Then \C 2 / R 2 \ ~ i]~ e . In the 
limit kr\ ^> 1, only the first term survives and \C 2 /R 2 \ oscillates at nearly constant amplitude. It is also interesting 
to note that the subhorizon perturbations evolve similarly to the Hamiltonian (|53|), ie. with constant amplitude. The 
superhorizon evolutions, on the other hand, differ from the Hamiltonian time dependence. Superhorizon perturbations 
are coupled to the backgound expansion and, in this case, the expansion is driven by density perturbations. So it is 
not surprising to find a scaling ~ r]~ 6 similar to that of decaying density perturbations. 

Finally we present results for the Bianchi type IX metric ( |70| ) , although due to the complexity of the Weyl tensor, 
we do not write out C 2 here. Because Bianchi IX is a vacuum solution with i? M „ = 0, we compute only |C 2 |, shown 
in Fig. H using the same initial data as in Fig. ||. For comparison, we also show f2 4 (introduced to bring out the 
structure at the scale of variations in |C 2 |). Notice that although C 2 oscillates (the kinks evident in |C 2 |, and which 
correlate with the peaks in f2, are points where C 2 becomes negative), the absolute magnitude diverges exponentially 
as the singularity is approached. The rate of divergence can be estimated from the "free fall" part of the trajectories 
during which a ~ (3 ~ e~ 3a , and the dominant terms in the square of the Weyl tensor scale as |C 2 | ~ e~ 12a for a«0. 

The above results continue to throw doubt on the utility of the C 2 definition of entropy. Indeed, the simple C 2 
measure seems to be again ruled out because of its inability to handle both the decaying and growing modes in a 
sensible fashion. 

Our results also point to important differences between the phase space and C 2 measures of entropy, as well as 
several other functions one might consider. As can be seen from above, C 2 is a local quantity, which will vary from 
point to point. As such it is not a useful measure of the global properties of spacetime, unless some sort of spatial 
average is introduced. By the same token, even though they are gauge invariant to first order, one can rule out the 
metric perturbations $ and h for the density and gravitational wave perturbations. These are also local quantities. 

The Hamiltonian in our examples could be considered on its own to be a measure of inhomogeneity since it is 
summed over the spatial coordinates and has a sensible time dependence. In regard to monotonicity in the time 
dependence, f2 appears to offer no advantages over H (except perhaps in the case of Bianchi IX, where we found that 
along the plateaus of constant H, f2 increased in the direction of increasing anisotropy). However, H alone does not 
provide a statistical description of a system in that it can be changed by the addition of an arbitrary phase, fi, on 
the other hand is a truly global quantity that expresses the entire allowable dynamical range equally for each of the 
oscillators in the spacelike hypersurfaces. It is not restricted to a particular phase realization, unlike any combination 
of variables constructed from metric components, which is. il thus presents the advantage over C 2 , H or any single 
solution to the differential equations in that it is global, shows a sensible time dependence and reduces to familiar 
entropy under appropriate circumstances. This is apparently true even in one area we have not yet addressed. 
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XI. CONNECTION WITH BLACK HOLES 



One of the questions one naturally wishes to answer is whether the entropy we have defined results in the well-known 
entropy of black holes. To establish the connection would strengthen any claim that the entropy function of this paper 
is in fact entropy. We now give a Bekenstein-style argument [^lj that the logarithm of the phase space does reduce 
to the entropy of black holes in the appropriate circumstance. The argument resembles the one we gave in 

§0 for 

the EM field and is also somewhat similar to one found in Zurek and Thorne |2^| ; we have, however, not seen this 
demonstration elsewhere. 

§[I] we showed that the phase space of harmonic oscillators, Eq. (|l3|), gives the classical limit for Einstein's 
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formula and results in a reasonable expression for the entropy of the electromagnetic field. This phase space was, 
with a slight change in notation, 

where lo is the angular frequency. 

Suppose we wish to construct a black hole out of photons, i.e., quantum oscillators. To do this, we must squeeze the 
oscillator system to within a Schwarzschild volume and the total energy of the oscillator system should equal M, the 
mass of the black hole. The latter condition implies that H = M = Ne, where e is the average energy of an oscillator. 
We also have e = v = 1/A, where A is the wavelength of the photon The minimum energy per oscillator needed 
to construct the hole corresponds to the longest allowed wavelength, which should be of order the Schwarzschild 
diameter, or A = AM. Let us, however, parameterize the wavelength as A = fM. Hence e = l/(/M) and 

N = fM 2 . (95) 

With the above expression for N , lo — 2ire = 2it/ fM and Stirling's formula, we quickly find by taking the logarithm 
of (@) that 

S = lnn = fM 2 . (96) 

The exact value of S therefore depends on /. A priori, we expect A ~ AM or / = 4. However, if Ap ~ 1/A, the 
uncertainty relationship implies that A may be as large as lQirM. In the former case we are a factor of 2tt 2 lower than 
the Bekenstein-Hawking |23) value of 87r 2 M 2 , in the latter case a factor of n/2 lower. Alternatively, if one chooses 
A = (2Ty 1 where the black-hole temperature T" 1 = 16tt 2 M = dS/dM ||] ,one recovers the exact result / = 8?r 2 . 

One might object that we have basically given a dimensional argument. Nevertheless, that the phase space of 
harmonic oscillators comes so close to the accepted result is striking. With hindsight, the phase-space approach 
makes clear that lnfi m N, so black-hole entropy must be of order M 2 . A Hamiltonian modified for quantum 
mechanical systems would, we expect, reproduce the usual result. Note, however, that unlike the cosmological models 
we have considered, the Hamiltonian here is not the ADM Hamiltonian for the black hole itself. The Hamiltonian for 
the Schwarzschild metric would presumably result in zero entropy since the canonical momenta are zero in the static 
case. Thus the harmonic oscillator Hamiltonian must be regarded either as perturbations on the background or as 
the Hamiltonian for the infalling oscillators; this latter corresponds to the usual approach for calculating black- hole 
entropy. We will explore these matters further and attempt a quantum mechanical calculation in a future paper. As 
it stands our current result shows that black-hole entropy can be treated profitably as a classical quantity. We also 
emphasize that, in contrast to the cosmological case, the black-hole Hamiltonian is easily interpreted as the energy 
and it is constant; the resulting phase space should then be the usual one. The main leap, evidently, in accepting the 
function we have termed gravitational entropy as genuine entropy lies not in the classical treatment, but in the use 
of time-dependent Hamiltonians. 



XII. FUTURE WORK 



We have evaluated the phase space for a number of models in the perturbative limit under the assumptions that: 1) 
the phases of the various components can be ignored; 2) that the system can be defined on spacelike hypersurfaces with 
some prescription for choosing boundaries; 3) the system is constrained by a Hamiltonian on each hypersurface. Under 
these assumptions In f2 appears to be a reasonable entropy function in that it increases with increasing inhomogeneity 
and not otherwise. Because the phase space for the perturbative spacetimes we have considered is computed using 
gauge invariant functions, entropy as we have defined it is thus also gauge invariant to first order. Moreover, it 



18 



can be identified with the entropy of more familiar situations. We also point out that the generalized second law 
of thermodynamics appears to be automatically satisfied. The generalized second law states that the sum of the 
thermodynamic and gravitational entropies in a closed, isolated system should always increase. Unless for some 
reason an increase in gravitational entropy actually causes a decrease in thermodynamic entropy the generalized 
second law should not be violated. (For some time this was not obvious in the black hole case, in which the hole can 
decrease the surrounding entropy (at the expense of increasing its surface area.) However, since our entropy becomes 
black hole entropy, this situation is evidently taken care of.) A more detailed investigation of this question may be 
warranted. 

We reiterate that all our calculations have been performed in the classical limit. We will present a quantum 
calculation for the black hole case in a future paper. 

Full, nonperturbative ADM calculations for inhomogeneous model systems would also be desirable. One system 
to examine is spherically symmetric collapse. However, in this case (as in classical orbital problems) the canonical 
coordinates and momenta appear to be coupled, making it difficult to perform the integrations. If the system is 
tractable, it may be possible to get black hole entropy by calculating the phase space available to a collapsing star or 
dust shells. 

These are a few problems we hope to examine in future work. The phase space approach is a generic one, applicable 
to a wide range of systems, including dust, radiation, N— body simulations, Newtonian and relativistic problems. 
Hence, the cases we have mentioned are probably only a small subset of those that can be examined. The more 
important message is that a consideration of the phase space available to general-relativistic systems appears to open 
a direct connection to statistical mechanics. This connection is well worth investigating. 



ACKNOWLEDGMENTS 



T.R. wishes to thank G.F.R. Ellis for his suggestion and great encouragement to work on the problem of gravitational 
entropy. He is also grateful to Richard Matzner and the members of the Center for Relativity for their hospitality 
during the early stages of this work. Thanks as well to Dilip Kondcpudi for a statistical mechanic's point of view. 
Both of us are grateful to Sherman Frankcl and the University of Pennsylvania physics department for their hospitality 
in the closing stages of this work. We thank Steven Brandt for his help in using Set Tensor, a script file he developed 
to interface with MathTensor and Mathcmatica. 



[1] Penrose, R. (1989), in Fourteenth Texas Symposium on Relativistic Astrophysics, ed. E. J. Fenyves (New York Academy 
of Sciences, New York). 

Wainwright, J. and Anderson, P.J. (1984), Gen. Rel. Grav. 16, 609. 
Goode, S.W. and Wainwright, J. (1985), Class. Quantum Grav. 2, 99. 
Goode, S.W., Coley, A.A. and Wainwright, J. (1992), Class. Quantum Grav. 9, 445. 
Bonnor, W.B. (1987) Physics Letts. A 122, 305. 
Smolin, L. (1985) Gen. Rel. Grav. 17, 417. 
Hu, B. and Kandrup, H. (1987) Phys. Rev. D 35, 1776. 

Brandenberger, R., Mukhanov, V. and Prokopec, T. (1993) Phys. Rev. D 48, 2443. 

Arnowitt, R., Deser, S. and Misner, C.W. (1962), Gravitation: An introduction to Current Research, ed. L. Witten, (John 
Wiley, New York) . 

Aarseth, S.J., Gott, JR. and Turner, EL. (1979) Ap. J. 236, 43. 

Efstathiou, G, Davis, M., Frenk, C.S. and White, S.D.M. (1985) Ap. J. Suppl. 57, 241. 
Anninos, P. and Norman, ML. (1996) To appear in Ap. J. 

Abramowitz, M. and Stegun, LA. (1972), Handbook of Mathematical Functions, (Dover Publications Inc., New York). 
Reichl, L.E., (1991) A Modern Course in Statistical Physics, (University of Texas Press, Austin). 
Reif, F. (1965), Fundamentals of Statistical and Thermal Physics, (McGraw-Hill Publishing Company, New York). 
Wald, R.M. (1984), General Relativity, (University of Chicago Press, Chicago). 
Mukhanov, V., Feldman, H. and Brandenberger, R. (1992) Phys. Rep. 215, 203. 

Misner, C.W., Thorne, K.S. and Wheeler, J. A. (1973), Gravitation, (W.H. Freeman and Company, New York). 
Ryan, M.P. and Shepley, L.C. (1975), Homogeneous Relativistic Cosmologies, (Princeton University Press, Princeton). 
Barrow, J.D. (1982) Phys. Rep. 85, 1. 
Bekenstein, J.D. (1973) Phys. Rev. D 7, 2333. 



19 



[22] Zurek, W.H. and Thorne, K.S. (1985), Phys. Rev. Letts. 54, 2171. 
[23] Hawking, S.W. (1976) Phys. Rev. D 13, 191. 



FIG. 1. Phase space trajectories for the two particle model, assuming one particle to be at rest (or equivalently more 
massive than the other) and mi = Grri2 = 1 with rt%2 mi. Three different trajectories are displayed: a reference curve of 
intermediate energy Eo and softening parameter ro, and two other curves varying Eo and ro independently to increase the 
phase space volume relative to the reference curve. 

FIG. 2. Contour plot of the Bianchi type IX potential V. Seven level surfaces are shown at equally spaced decades ranging 
from 1CP 1 to 10 . For V > 1, the potential is open and exhibits a strong triangular symmetry with three narrow channels 
extending to spatial infinity. For V < 1, the potential closes and is approximately circular. 

FIG. 3. The Hamiltonian and entropy for Bianchi type IX as a function of a. The evolution is initialized at a — with 
the following data: f3+ = (3- = 0, $+ = 2, and /3_ = 1, and run both forward in time and backward towards the singularity 
a — ¥ — oo. 

FIG. 4. A comparison plot of the phase space volume (represented by Q 4 ) and the magnitude of the Weyl tensor squared 
\C a p S C al ^ s \. The kinks evident in the |C 2 | curve represent regions where C 2 < 0, which are correlated with the regions in 
which Q drops after reaching a peak value, ie. the wall collisions. 
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